MU'KiVOH  msOmnON  H SI  I'HAKI 

NAIK’NM  Ul'NlAl'  y'l  UNt*AKl'S  \ 


ADA072153 


? 

\ 

i 


THE  APPLICATION  OF  THERMAL  AND  CREEP  EFFECTS  TO 
THE  COMBINED  ISOTROPIC-KINEMATIC  HARDENING  MODEL 
FOR  INELASTIC  STRUCTURAL  ANALYSIS  BY 
THE  FINITE  ELEMENT  METHOD 


LEVEL1 

aerospace 


(ffr) 


< 


n 


enqineeRinq 

department 


TEXAS  A&M  UNIVERSITY 


David  H.  Allen  and  Walter  E.  Haisler 


Iff  ice  of  Naval  Research 
Department  of  the  Navy 
Washington,  D.  C.  20025 


PERIMENT  ST 

- '■  . , m. 


THE  APPLICATION  OF  JHERMAL  AND  £REEP  EFFECTS  TO 
THE  COMBINED  pOTROPIC-KINEMATIC  HARDENING  JNODEly 
FOR  JNELASTIC  JTRUCTURAL  ANALYSIS  BY 
THE  FINITE  ELEMENT  METHOD  „ 


David  H./Allen  art  Walter  E./Halsler 


Aerospace  Engineering  Department 
Texas  A&M  University 
College  Station,  Texas  77843 


Office  of  Naval  Research 
Department  of  the  Navy 
Washington,  D.  C.  20025 


Contract  No/  Nflfl 
Task  No.  NR 


Technical  Report  Nd\' 3275-79-3 


W i 


V 


THE  APPLICATION  OF  THERMAL  AND  CREEP  EFFECTS  TO 
THE  COMBINED  ISOTROPIC-KINEMATIC  HARDENING  MODEL 
FOR  INELASTIC  STRUCTURAL  ANALYSIS  BY 
THE  FINITE  ELEMENT  METHOD 

David  H.  Allen*  and  Walter  E.  Haisler** 
ABSTRACT 


Accession  For 


NTIS  GfiA&I 
K>C  TAB 
Unsaa junced 

Justification 


By_ 

"Ibuttry/ 

Avnllabinty^codf 

Avail  and/or 
special 


Dist 

ft 


In  this  paper,  a formulation  is  presented  for  a combined  isotropic- 
kinematic  hardening  rule  which  includes  both  creep  and  thermal  effects. 

The  resulting  model  accounts  for  temperature  dependent  material  properties 
as  well  as  the  Bauschinger  effect  for  reverse  loading,  ^^nrtual  work 
equation  and  resulting  matrix  equations  are  derived  for  implementation 
of  the  model  into  a general  purpose  finite  element  code^ 

INTRODUCTION 


In  the  past  decade  considerable  research  requiring  structural  analysis 
including  thermal  and  creep  effects  has  been  done  in  fields  such  as  laser 
technology  and  nuclear  reactor  analysis.  Thus,  it  is  desirable  to  modify 
existing  nonlinear  structural  analysis  codes  to  encompass  thermal  and  creep 
effects . 

Researchers  at  the  Massachusetts  Institute  of  Technology1  as  well  as  the 

2 

Oak  Ridge  National  Laboratory  have  included  thermal  and  creep  effects  in 
structural  analysis  programs  for  the  cases  of  isotropic  hardening  and  kine- 
matic hardening  applied  to  metal  plasticity.  However,  as  has  been  pointed 
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out  by  Haisler  , there  may  be  certain  material  types  and  stress  histories 
for  which  these  models  are  inadequate. 

One  model  which  has  been  shown  to  predict  cyclic  loadings  quite  well 

for  metals  under  Isothermal  conditions  is  the  combined  isotropic-kinematic 

4 

hardening  rule  . This  paper  outlines  the  formulation  used  to  obtain  the 

finite  element  equations  necessary  to  Include  thermal  and  creep  effects  in 

a combined  isotropic-kinematic  hardening  model  currently  in  the  general 

5 

purpose  code  AGGIE  I . 

In  the  first  section  of  this  paper,  the  combined  isotropic-kinematic 
workhardening  rule  will  be  utilized  to  obtain  a material  model  which  accounts 
for  creep  and  thermal  effects.  Next,  the  derivation  of  the  linear  virtual 
work  equation  within  a total  Lagrangian  description  is  presented  for  the 
material  model  described  above.  The  constitutive  law  will  then  be  applied 
to  the  virtual  work  expression  and  the  resulting  relation  will  be  discretized 
to  obtain  a matrix  formulation  which  is  applicable  to  finite  element  modeling. 
The  iterative  technique  necessary  to  obtain  equilibrium  for  the  nonlinear 
virtual  work  equation  will  then  be  presented.  Finally,  the  constitutive 
law  will  be  specialized  for  certain  load  conditions  and  the  resulting  matrix 
equations  will  be  given. 

THE  CONSTITUTIVE  RELATION 

In  order  to  determine  the  quantities  in  the  constitutive  law,  it  is 
necessary  to  employ  an  appropriate  material  modeling  theory.  The  model 
used  In  this  report  is  the  incremental  theory  of  plasticity.  In  this 
method  a yield  criterion  such  as  the  Von  Mises  yield  condition  is  used 
with  an  associated  flow  rule  and  a so-called  workhardening  rule. 

The  workhardening  rule  utilized  here  is  the  combined  isotropic- 
kinematic  hardening  rule,  which  accounts  for  both  expansion  and  translation 
of  the  yield  surface  in  stress  space  and  is  described  by  the  yield  function 
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where 


(i) 


F<su  - «ij>  ■ k‘*'<r  T» 

j » coordinates  of  the  yield  surface  center  In  stress  space, 

<}  - plastic  strain  Increment  tensor, 

T ■ temperature, 

and  /dE^j  represents  the  plastic  strain  history  dependence  In  the  yield 
function. 

The  absence  of  temperature  dependence  In  on  the  left-hand  side  of 

the  above  statement  signifies  that  there  Is  assumed  to  be  no  translation 

of  the  yield  surface  due  to  temperature  changes.  Hence,  thermally  Induced 

yield  surface  modifications  are  assumed  to  be  strictly  Isotropic.  Phillips10 
has  shown  In  his  experimental  work  Involving  metals  that  although  thermally 

Induced  yield  surface  modifications  may  be  slightly  anisotropic  for  highly 

distorted  yield  surfaces,  the  absence  of  translation  Is  verified  for  any 

yield  surface  In  stress  space.  An  Identical  assumption  Is  used  In  the 

kinematic  hardening  rules  employed  at  the  Oak  Ridge  National  Laboratory 

and  the  Massachussets  Institute  of  Technology. 

From  the  above  yield  function  It  can  be  seen  that  a statement  of 

consistency  during  plastic  loading  Is  given  by 

&1J  dS«  ‘ d“'J  * *“  ('<IE'J)  ■ dT  ■ 0 (2, 

where  dk(/dE^)1s  defined  to  be  the  change  In  the  yield  surface  size  with 
respect  to  the  history  of  plastic  strain.  Because  the  above  consistency 
condition  Is  difficult  to  employ  directly,  the  consistent  statement  employed 
by  Ziegler  for  Isothermal  kinematic  hardening  Is  modified  and  utilized  here. 
The  applicability  of  the  Ziegler  modification  to  combined  Isotropic-kinematic 
hardening  has  been  previously  established  by  this  author7.  As  employed  In 
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the  combined  Isotropic-kinematic  hardening  rule  the  Ziegler  modification  Is 


fiy  ‘<J  ■ Sy  ^Ij  + ^ (^Etj)  (J) 

where  c Is  a scalar  constant  (called  the  hardening  modulus)  to  be  determined. 
Substituting  the  above  relation  into  equation  (2)  gives 

K dS” ' K cdt'J ' dT  ■ 0 (4) 

The  above  consistency  condition  has  the  added  advantage  that  It  cap  be 

t 

Interpreted  as  a statement  controlling  nonlsothermal  loading  and  unloading. 
Since  during  neutral  loading  the  plastic  strain  Increment  Is  zero,  a statement 
representing  plastic  loading  would  be 

f£  dSn  - 2kUdT  >0 

3S1j  1j  9T  (5) 

and  It  follows  that  a statement  governing  unloading  Is 

fftJ  dS1J  - 2kff  <"  <°  («) 

It  can  be  seen  that  thermally  Induced  changes  In  the  yield  surface  size 
can  be  significant  In  affecting  the  loading  process.  For  example,  equation 
(5)  correctly  predicts  loading  when  the  stress  Increment  Is  zero  and  the 
temperature  Increment  Is  positive  since  3k/aT  Is  negative  for  most  metals. 

In  fact,  a small  negative  stress  Increment  may  even  produce  loading  If  the 
temperature  Increment  produces  a large  enough  contraction  of  the  yield  surface 
to  ensure  that  equation  (5)  remains  positive. 

The  total  stress  at  time  t for  a material  with  no  strain  rate  dependence  Is 


where 


* DE  |r^  _ pPt  _ Fct  pTt\ 

b1j  D1jmn  'Sin  Emn  Emn  “ Emn> 
l 

°1jmn  " el*st^c  modulus  tensor, 

Ct 

■ creep  strain  tensor,  and 


« thermal  strain  tensor. 
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The  total  differential  of  the  above  equation  Is 

dSU  * °U  * d0tJ»n<C-C-C-C>  »> 

where  the  Increment  dDjjmn  Is  caused  by  a temperature  change. 

The  above  equation  Is  exact  only  In  the  limit  as  the  stress  Increment 
approaches  zero.  In  spite  of  this  fact,  many  researchers  have  used  equation 
(8)  to  define  the  Incremental  constitutive  relation  for  finite  load  and 
temperature  steps.  Unless  the  material  properties  vary  significantly 
with  temperature,  the  above  approximation  Is  satisfactory,  but  for  the  case 
where  material  properties  depend  heavily  on  temperature  (e.g.,  above  half  the 
material  melting  point  In  many  metals)  an  additional  term  should  be  added 
to  account  for  finite  step  sizes.  To  obtain  this  term,  consider  the  stress 
at  time  t+At,  which  Is  given  by 

-t+At  nt+At  /pt+At  pPt+At  pCt+At  pTt+Atv 

51j  u1jmn  ltmn  " fcmn  “ bmn  " bmn  ' (9) 

The  change  In  stress  during  a load  step  Is  then  given  by  the  difference  In 
equations  (9)  and  (7): 


Dt+At  fct+At  _ pPt+At  pCt+At  pTt+Atx 
u1jmn  'Lmn  " Snn  Sin  " Snn  ' 


- D1Jmn  <E1  * C - C - O 
Equation  (10)  may  be  written  equivalently  as  follows: 


. dS1j 


(D 


:P 

‘mn 


<2> 


'Ijmn  +‘  d01Jmn^  ^Emn  + dEmn  ~ Emn  * dEn 

■ C - dEL  - - «L>  - »!*■  <C  - C - *2 

A further  refinement  of  equation  (11)  yields 

dS1J  ■ D1Jmn  <dE»n  * dC  * <„  - dO 
* dD1jmo  <E1  - C - E*J  - O * dD1Jmn  <dE*n-dEL-di-dEL.) 


(10) 


(11) 


(12) 
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A comparison  of  equations  (8)  and  (12)  reveals  that  the  differential  form 
neglects  the  second  order  effect  produced  by  the  last  term  In  equation  (12). 
Our  experience  shows  that  the  second  order  term  may  be  neglected  In  all 
cases  except  one;  namely,  that  the  temperature  change  during  a load  step 
produces  an  elastic  constitutive  tensor  modification  which  will  cause  the 
last  term  In  equation  (12)  to  become  significant  when  compared  to  the  first 
two  terms.  Since  we  have  already  encountered  test  cases  where  this  term  Is 
slgnlflcant.we  have  Included  It  In  our  formulation.  We  have  chosen  for  pro- 
gramming purposes,  however,  to  write  equation  (12)  in  the  following  form: 

dS, 


(dEmn  ' dEL  - dEL  - dET  ) 
1 Jmn  mn  mn  mn  mn' 


rPt 


ECt  _ ETtj 
mn  mn' 


+ dD1jim/Emn  “mn 
Now  recall  that  the  normality  condition  requires  that 


(13) 


dE 


Ij 


dx 


3F 

laS 


ij 


(14) 


Substituting  this  relation  into  equations  (4)  and  (13)  gives 


9F 

55. 


dS 


aF 

ij  " 3$, 


cdx 


and 


' 2kff  dT  =0 


’Ij  M1j 

dS1j  " Dfjmn  ^dEmn  " dX3S  " dEmn  “ dEmn^ 

mn 


05) 


mn 

* "tjmn  - C ' 


O 


(16) 


Next,  substitute  equation  (16)  into  equation  (15)  and  solve  for  dx: 

dx  * nt+At  dE  - iL  n^+At  jpC  _ JF  nt+At 

L*ij  ijmn  mn  a5^j  Ijmn  mn  as^j  ijmn  Emn 

+ dD  _ 3F  jn  pPt  aF 

as,,  ijmn  Sin  a5..  auijmn  Sin  " at..  dD..  Ect 
■J  ij  ijmn  mn 

££  dD, , ETt  - 2k^dll /\c  Ur  H + nt+At  if.  & 1 
3S1j  Ijmn  mn  at  TJ/  L at^  a?^.  Dijmn  aS^  asmn J 


(17) 
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Now  substitute  equation  (17)  back  Into  equation  (16): 


( if.  Qt+At  dE 

1 jmn  l atmn  3?^  utuvw  atvw 


3F  nt+At  .rC 
3$*,.  utuvw  atvw 

vil 


jt+At 

;tuvw 


+ 


dO 


tuvw 


- •?£.  Hn  cpt  _ 3_F  dn  cCt  jr*  _Tt 

3Stu  dDtuvw  Evw  3Ttu  dDtuvw  Evw  " 3Stu  dDtuvw  Evw 

- 7V^L  /(r  3F  3F  + nt+At  3F  3F  \ .cC  .CT  1 

2k5T  «)/{c  5?pq  5?pq  “pqrs  5Jpq  5 ?J  ' dEmn  * dEmn  J 


+ dD4<  (E*  - FPt  - ECt  - E™) 
Ijmn  v mn  'mn  cmn  mn; 

Rewrite  equation  (18 ) In  the  following  form 


(18) 


dS 


ij 


_t+At  3F  3F  Dt+At 

^t+At  ^vw  ^vw  ^tu  tLimn 
51  jmn  ‘ ~TF  a?  . "t+At  3F  a? 

^nn  PR'S  3^pq  3S 


pr  pq 

dD* 


t+At  aF  3F 

°tjvw  gv»  *» 

TaF  aF  . "t+At  3F  aF 
V ^pq^pq  pqrs^pq^. 


rs 


rs, 


- E, 


(dE, 


Pt 

mn 


mn 


- dEL  ~ dE*  ) 

mn  mn 


ECt  _ ETtj 
mn  mn ' 


nt+At  9.  3k  3F 
°1  jmn  2k?f  3S„ 


mn 


TaF  aF  TTf+T t aF — aF 
^ Vs  Spq  Srs 


dT 


or 


+ dO^  (E*  - EPt  - ECt  - ETt  ) 
ijmn  vmn  mn  cmn  mn 


dSiJ  ■ C5jmn  <dEmr,  ’ dEmn  ' dO  + dP1j 


09) 


(20) 


where,  noting  that  the  elastic-plastic  constitutive  tensor  Is  superscripted 
t because  It  Is  obtained  by  using  the  hardening  modulus  at  time  t. 
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,t*M  3F  JF  „t«t 


C*,  - 0?*“  - (J,“  Hvw  *tU  tum" 

1ji.fi  Umn  ";T'~'aT  .“tHtaF  aF 

C ^nn  ^nn  pqrS  35__  ?? 

pq  pq  ^ pq  rs 


t +At  aF  aF 
Uvw  3?yw  3?t(j 


l3r  aF  , nt+AraF  aF 

c ainft  f5_„  upqrs  3?  5!f 
pq  pq  pq  rs 


- C ■ 4n  - O 


nt+At  ,.ak  aF 
jmn  2kaT  a? 


71F  aF  . nt+At  aF  aF“ 

C ^rS 

pq  pq  ^ pq  rs 


dT  + dD4 . (Et  - EPt  - ECt 


Ijmn  ' mn 


mn  mn ' 


It  Is  of  particular  Interest  to  note  that  In  the  case  of  isothermal 


loading,  dP^  is  zero  and  the  model  derived  by  Hunsaker  for  Isothermal 
combined  isotropic-kinematic  hardening  is  recovered. 


The  last  term  in  equation  (20)  arises  solely  due  to  the  admittance 


of  nonisothermal  action.  The  terms  in  equation  (22)  which  contain  dDtumn 


represent  the  stress  Increment  due  to  a thermally  Induced  change  in  the 


elastic  constitutive  tensor.  The  term  containing  3k/3T  in  the  same  equation 


accounts  for  the  stress  increment  caused  by  a modification  in  the  yield  surface 


due  to  a temperature  change. 


In  order  to  better  understand  the  term  dP^  in  the  constitutive  relation. 


consider  the  case  where  a member  is  loaded  such  that  the  total  elastic  strain 


tensor  is  non-zero  and  the  total  plastic  strain  tensor  is  zero.  Now,  if  a 


load  Increment  accompanied  by  a temperature  change  is  Introduced  such  that 


the  elastic  strain  increment  is  zero,  a stress  increment  will  occur  which  is 


equivalent  to  the  product  of  trte  change  in  the  elastic  constitutive  tensor 


and  the  total  elastic  strain  tensor. 


In  order  to  completely  define  the  above  constitutive  relation  it  is 
necessary  to  account  for  sk/3T,  c,  and  u.  To  obtain  the  first  quantity 


- 1 
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note  that  3k/aT  Is  the  yield  surface  gradient  with  respect  to  temperature. 
This  value  may  be  obtained  by  a suitable  Interpolation  of  Input  uniaxial 


stress-strain  data. 


To  evaluate  c,  first  solve  equation  (15)  for  c.  The  result  Is 


p dSj  • 

‘Ij  ij 
If — aT 


2klr  dT 


al__  aSnn 

pq  pq 


"asnn  5? 
pq  pq 


which  may  be  rewritten 


3F  HQ 

dS’J 


*a?T  as 

pq  pq 


2klrdT 


dsu 


Now  define 


3Snn 

pq  pq 


so  that  equation  (24)  may  be  written 


c-2h7,  2k^dT 
3 H I 1 ' sF  » 

\ *lj  *« 


Sinc,e  H'  Is  a stress  Invariant  it  may  be  determined  from  uniaxial  stress- 


strain  data  . 


The  last  term  in  the  above  equation  causes  some  computational  difficulty 


since  It  involves  the  stress  increment  tensor  which  is  itself  a function  of 


the  hardening  modulus  c.  This  shortcoming  has  been  circumvented  by  some 


researchers  by  simply  neglecting  it.  This  assumption  is  a clear  violation 
of  the  consistency  condition  and  can  cause  departures  of  the  stress  tensor 


from  the  yield  surface  during  loading.  Therefore,  the  thermal  dependence 


of  the  hardening  modulus  Is  not  neglected  herein. 
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In  order  to  satisfy  the  condition  that  c must  always  be  a positive 
quantity  for  a workhardening  material  (i.e.,  H'  Is  strictly  positive). 

It  Is  necessary  to  show  that  the  last  term  In  equation  (26)  is  always  less 
than  unity.  This  may  be  verified  quite  simply  by  rearranging  the  loading 
condition  given  In  equation  (5). 

Finally,  to  account  for  u,  recall  Ziegler's  rule, 

do1j  = *,(sij  - a1j)  (27) 

Substitute  equation  (27)  Into  equation  (2)  and  solve  for  p. 


««  - 2klrdT  - 2kdk  <'<> 

-1 it 


(s, 


mn 


otmn)  a5_ 


'mn  (28) 

The  thermal  term  in  the  above  equation  can  be  seen  to  be  the  correction 
in  the  yield  surface  center  translation  tensor  which  accounts  for  isotropic 
expansion  due  to  a temperature  Increment. 

This  completes  the  formulation  of  the  constitutive  law. 


THE  INCREMENTAL  VIRTUAL  WORK  EQUATION 


In  order  to  formulate  a statement  of  equilibrium  which  will  include 
thermal  effects,  first  consider  the  virtual  work  expression  within  a total 
Lagrangian  description  as  previously  described  by  Hunsaker®: 

/sfr  5Eijt+4t  dv  ■ /TkVA  *AFk% dv 

Vo  A Vo  (2! 

where 

-t+At 

a1j  * the  2nd  Plola-Kirchhoff  stress  tensor  at  time  t+At  referred 
to  the  Initial  configuration  VQ  at  time  t = 0, 

6E^At  = the  variation  In  the  Green-Lagrange  strains  at  t+At  referred 
to  the  Initial  configuration  VQ, 
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the  surface  tractions  at  time  t+At  referred  to  the  surface 


of  the  configuration  A, 

<3uk  ■ the  variation  In  the  displacements, 

p0  * the  local  density  In  the  Initial  configuration,  and 

Fk  * the  body  force  per  unit  mass  at  time  t+At  referred  to 
the  Initial  configuration  VQ. 

Since  the  Inclusion  of  thermal  effects  Into  the  system  will  affect  only 
the  left  hand  side  of  equation  (29),  make  the  following  simplifying  definition: 


.t+At  . 


/ VukdA  + / °oFk% 


As  oucllned  by  Hunsaker,  In  order  to  formulate  an  Incremental  stiffness 


procedure,  first  let 


S1 3 + AS1j 


pt+At  _ rt  , .r 

E1j  " Eij  + AE1j  ( 

where  and  are  respectively  the  2nd  Plola-KIrchhoff  stress  tensor 
and  Green-Lagrange  strain  tensor  at  time  t referred  to  the  Initial  config- 
uration, and  aS.j j and  AE^j  are  stress  and  strain  Increments. 

The  virtual  work  expression  [equations  (29)  & (30)]  thus  becomes 


y<su + 4su> { <E?j + 4Eu> dv  ■ 


Since  the  strains  at  time  t are  known  they  need  not  be  varied  and 
equation  ( 33)  reduces  to 


Mi  * 4sij> s <4Eij> dv 


] 


Now  decompose  the  strain  Increment  into  linear  and  nonlinear  components: 


where 


iEiJ  ■ 

lEij 

+ AFNL 

+ AE1j 

(35) 

3AUj 

3uk* 

3AUk  SAU^  au^X 

9Xi 

3X1 

*\T  + 33TT^~/ 

(36) 

NL  _ 1 / 

9Aui 

3AUk\ 

■1J  "?\ 

3X.. 

3Xj  ) 

(37) 

and  Au.j  are  the  displacement  Increments  from  time  t to  t+At.  Equation 
(34)  becomes 

f (Sjj  + AS^)  6 (AEjj  + AE^)  dV  = «Wt+At 
Vo  (38) 

Next  assume  the  stress  Increment  during  any  time  step  may  be  written 

aC  = />r  _ *rC  ,rT  \ . »n 


AS1j  “ Cijmn  ^AEmn  ' AEmn  ' AEmn^  + APi 


where 


Cijmn  = tfie  effectlve  modulus  tensor  at  time  t, 
r 

AEmn  * the  creep  strain  increment  from  time  t to  t+At, 

AEinn  = the  thermal  strain  increment  from  time  t to  t+At,  and 

AP^j  ■ the  stress  increment  due  to  creep  and  thermal  loading. 

It  can  be  seen  that  the  above  constitutive  relation  is  consistent 
with  the  material  law  presented  In  equation  (29). 

Now  substitute  equation  (35)  into  (39)  and  (39) into  equation  (38) 
to  obtain 
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/su  '4E<J  dV  + /SU  54E!S  dV  4 fAm  < “E(j  * 

Vo  ».  V0 

+ /CU  <„  “E7j  dv  ■ f^Sm  A<ln  “E7j  * 

Vo 

“EU  dV'/C‘ J«n‘4,“E1JdV 
Vo  vo 

-/C?jm.  4EL  “E1j  dv  -/«**»  4E1  “eij  «» 

Vo  »o 

-/ C1  j™  aEL  S4E1  j dV  ♦ /“’ll  54Eij  «» 


1 

I 


6AE^  dV  « aWt+At 


0 (40) 

Terms  4,5  and  6 In  equation  (40)  are  nonlinear  In  the  displacement 

Increment.  In  order  to  allow  for  a linear  Incremental  formulation,  It  will 
be  assumed  that  the  displacement  Increment  Is  small  so  that  these  terms  can 
be  neglected  at  this  time.  However,  this  approximation  will  be  accounted  for 
later  In  the  section  describing  equilibrium  Iteration.  With  this  modification, 
the  virtual  work  expression,  with  equation  (30)  Included,  simplifies  to 
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J° oV“k  dvt/c?jn,n  4EL  “Eij  dv  */Ai  sae75  dv 


-/c1j«n  a4  “Eu  dV  »4  S4E"i  dv 


*fru  “e1j  dV  -f\ “k"-/s‘j  “Eij dv  +/CU 4E™  “Etj dV 


A 


*/ C1jmn  4E™  “E1j  dv  - /4P1J  S4E1.j  dv 


o vo  (41) 

The  above  virtual  work  equation  applies  for  elastic-plastic  action 
as  well  as  creep  and  thermal  effects.  Note  that  If  thermal  and  creep  effects 
are  neglected  equation  (41)  reduces  to  the  virtual  work  expression  presented 
by  Hunsaker. 

FINITE  ELEMENT  DISCRETIZATION 

The  process  of  obtaining  a finite  element  approximation  of  the  linearized 
virtual  work  expression  [equation  (41)]  will  now  be  summarized.  For  a more 
detailed  description,  the  reader  is  referred  to  Hunsaker. 

First  recall  that  for  an  Isoparametric  element  the  displacements,  u^ , 

and  coordinates,  x^,  of  points  In  a body  may  be  related  via  the  element  nodal 
displacements,  u“,  and  coordinates,  x“,  by  using  shape  functions,  h“. 


.a  a 

ui  = h U| 


(42) 


.a  a 
xi  ■ h Ui 


(43) 


LI 
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Note  that  repeated  Greek  Indices  Indicate  summation  over  the  total 
number  of  nodes  In  the  element. 

The  strain  Increment  may  thus  be  related  to  the  displacement  by 

4Eij  ■ 7 + s'  4UJ  + C fxj  “k  < * 4) 


(44) 


AFNL  . 1 3h°  3hP  . a 8 
AE1  j ■ 7 ax,  Jx.  Auk  Auk 

The  displacement  and  strain  variations  In  equation  (41)  may  now  be 
written  for  a single  element  as 


(45) 


6AE, 


6AE 


rs 


NL 

,J  9uk 


3AEE 

3AEE 

= — — 6U> 
5U“  1 

rs 

3Au“ 

3AEH1: 

6u, 


6uk  = 


9AU, 


(46) 


(47) 


3U 


6U. 


9U 


li  <SU* 

o 1 


_ 9Aut  „ 
° ~ i 6u“ 

a l 


9AU 


1 


Now  define  {Rt+At}  by 


(48) 


>r  -f\  ^ dAt/p°Fk2 

J 3AU?  J 3AU1 

A 1 t. 


dV 


(49) 


Next,  In  order  to  evaluate  the  second  Integral  In  equation  (41),  first 
define 

(AS)  = [C]  {AE}  (50) 

where 

(AS)T  * I AS 1 x AS22  AS33  ASj2  AS13  AS23 I (51) 
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and 


{AE}  ■ |_AE11  aE22  AE33  AE12  AE13  AE23J  (52) 

and  [ C ] Is  the  matrix  representation  of  Cijmn  written  In  engineering 
strain  notation  . 

Also,  write 


‘ tBL]  <4“«>  (53) 

where  [ BL  ] Is  the  linear  strain  displacement  transformation  matrix 
obtained  from  equation  (44).  Thus,  when  applied  at  all  nodes  In  the  element, 
the  second  term  In  the  virtual  work  expression  may  be  written 


1 3AE„_ 

cH„  ~ 

1JrS  3U? 


dV  - [KJ  {Aua} 


where  n Is  the  number  of  nodes  In  the  element  and 


(54) 


[Ki]  * / Cbl]  T [C] 


d V 


(55) 


The  remaining  terms  on  the  left  hand  side  of  equation  (41)  are  found 
In  a similar  manner  to  yield 
3 n 


k»l  a*l  I V k V 

•-ft  F A 


3AU, 


• ftjm  < « + 

V.  k 


•NL 


■/*»**] 


[Knl]  {AUa> 


(56) 
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where 


n W ' /%  ] T ([S]  - [C]  [4EC]  - [C]  [4ET]  + [4P])  [Bnl]  dV 

vo 

and  [Bnl]  Is  the  nonlinear  strain-displacement  transformation  matrix 
obtained  from  equation  (45). 

The  remaining  four  terms  In  equation  (41)  may  be  written 


(57) 


3 n 


k*l  a=l 


r t 3A^H  f t C 3A^11 

/s^  dvt/c^AE™^dV 


LV„ 


' 94U?  i 1J  S4U? 


dV 


j [Bl]  T({S>  - [C]  {AEC}  - [C]  {AET}  + {AP})  dV  *-{Ft} 


(58) 


If  the  body  force  term  In  equation  (41)  Includes  Inertial  effects, 
the  following  term  may  be  separated  out  and  brought  to  the  left  hand  side: 

t+At  3AU. 


/• 


pouk 


3AU 


£ dV  » [M]  {uot+At} 
1 


■0  (59) 

where  [M]  Is  the  consistent  mass  matrix. 

The  flobal  equations  of  motion  for  the  entire  structure  may  thus  be 
assembled  from  the  element  equations  to  obtain 


W {i0t4At>  ♦ ([•£)  ♦ [K^])  C4ua)  - (Rt+At)  - (Ft> 


(60) 


The  solution  of  the  above  equations  will  result  In  approximate  displace- 
ment increments  due  to  the  neglectlon  of  nonlinear  terms  In  the  virtual 
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work  expressions.  This  approximation  may  be  accounted  for  by  performing 
equilibrium  Iteration. 

EQUILIBRIUM  ITERATION 

Since  equation  (41)  is  not  exact,  It  Is  necessary  to  perform  equilibrium 
Iteration  on  each  load  step  to  obtain  the  exact  displacements.  This  Is 
accomplished  by  utilizing  the  virtual  work  equation  In  Incremental  form 
[equation  (33)]  and  recognizing  that  for  exact  equilibrium 


where  (R  ) Is  the  Internal  force  term  defined  by  equation  (49),  and  {F 
Is  the  external  force  term  defined  by 


For  dynamic  problems,  the  Inertia  effect  may  be  Included  In  equation  (62) 
To  Improve  the  Initial  estimate  one  can  define 


so  that  for  exact  equilibrium  (f)  = {0}.  Applying  the  modified  Newton  method 
to  solve  the  nonlinear  system  of  equations  given  by  (61),  one  obtains 


Notice  that  (F  Is  defined  by  equation  (62)  and  not  by  equation 
(58)*  It  should  be  noted  that  a strict  application  of  the  Newton  method  would 
require  the  use  of  [KLt+At]  and  [K**At].  However,  the  use  of  [kJ]  and  [kJl] 
generally  Is  sufficient  and  provides  a computationally  faster  solution.  If 
convergence  Is  slow,  the  stiffness  terms  may  be  updated,  l.e.,  use  [K*+At] 

rift+Atn 


SPECIAL  LOADING  CONDITIONS 


} 


The  equations  of  motion  described  by  equation  (60)  represent  an  elastic- 
plastic-creep-thermal  load  response  for  a material  with  temperature  dependent 
material  properties.  Under  certain  special  loading  conditions  the  constitutive 
law  and  the  equations  of  motion  may  be  simplified  so  that  the  solution  scheme 
may  be  computationally  Improved.  The  special  cases  considered  herein  are 
quasi -Isothermal  elastic  action,  nonlsothermal  elastic  action,  and  quasl- 
f so thermal  elastic-plastic  action,  where  quasi -Isothermal  action  is  defined 
to  be  any  load  response  which  Includes  thermal  loads  but  for  which  the  material 
properties  may  be  considered  constant  over  the  temperature  range  encountered 
during  the  load  history.  The  three  loading  conditions  above,  together 
with  the  previously  derived  nonisothermal  elastic-plastic  action,  are  chosen 
because  of  their  diverse  programmabilities  Into  the  computer  code  AGGIE  I. 

Quasi -Isothermal  Elastic  Action 


During  quasi -Isothermal  elastic  action  the  constitutive  relation  given 
by  equation  (39)  may  be  simplified  to 


AS,,  * (AEmn  - AE*L  - AE  T) 

1j  ijmn  ' mn  mn  mn  ' 

where  Is  the  elastic  constitutive  tensor  at  time  t. 

linearized  virtual  work  expression  Is 


For  this  case  the 


/Vk  % dV  ♦ Msm  4Eim  54E1j  dV 


/ ‘isij 


DU„n  - »1>„  O dV 


V“k  « 


- / <S?j  - “{j™  *4,  * °Tjmn  4E»T>  54EU  « 


If  certain  terms  are  defined  In  accordance  with  the  above  equation, 
the  equations  of  motion  given  In  equation  ( 60)  need  not  be  altered.  For 
the  loading  case  specified  here  the  newly  defined  terms  are 

1*^)  " / £bnl]  T - W Uec]  - COD  [aet])  [bnl]  dV 

V°  (68) 

and 

(F*)  « J [Bl]  T ({S>  - [D]  {aEC}  - [D]  {AET})  dV 

V°  (69) 

where  [D]  Is  the  elastic  constitutive  matrix.  In  addition,  [S]  and  {S}  are 

the  stress  matrix  and  stress  vector,  respectively,  described  by  Hunsaker, 

and  [E]  and  {E>  are  similar  strain  Increment  representations. 

Nonlsothermal  Elastic  Action 

When  the  load  response  Is  nonisothermal  elastic  the  stress  at  time 

t Is 


B1jmn  ^Emn 


- E_ 


Ct 


- E. 


Ttx 


1j  "ijmn  '~mn  ~mn  “mn 
The  total  Increment  from  equation  (13)  Is 


(70) 


AS1j  “ DlSn  (AEmn  ‘ ’ O 

+ AD1jmn  ^n  " Emn  " Emn  ^ (71) 

where  AD^j^  represents  the  change  In  the  elastic  constitutive  tensor  during 
a time  step  due  to  a temperature  change.  Note  that  the  entire  last  term 
represents  AP^j  In  equation  (41).  Since  D^mn  can  be  determined  directly 
from  Input  stress-strain  data,  the  above  constitutive  relation  may  be 
substituted  Into  the  virtual  work  expression  which  becomes 


-20- 


- / »oFk4uk  « * / Cn  aEin  **1J  « 


/re*  - aF  E nt+A^  af  t 

^Ij  D1jmn  AEmn  " °ijmn  AEmn 


Ct 


AB1jmn  ^Emn  “ Emn 
ro 

- C*n  jiEij « ■ A% dA  - / ts5j  - ®«2  4E„nc  • C; 

A Vo 

ADijmn  ^Emn  Emn  Emn  ^ 4aEij  dv  (72) 

Once  again  the  discretization  process  yields  equation  (60)  where  now 

B*]  * / Cbnl]  t ([s]  - [o]  tAE°D  -[o]  [aet] 


and 


+ [AO]  [E*  - ECt  - ETt])  [Bnl]  dV 


fpt}  55  / Cbl]  T ({S}  ' [D]  ueC}  * ^D]  (AET) 


(73) 


+ [AD]  (E*  - ECt  - LTt»  dV  (74) 

Quasi -Isothermal  Elastic-Plastic  Action 

When  the  load  response  Is  such  that  the  thermal  effects  may  be  considered 
quasi -1 sothermal , the  elastic-plastic  constitutive  relation  reduces  to 


ASij  Cijmn  (AEmn  “ AEmnC  ‘ AEmnT> 

With  this  simplification  the  virtual  work  equation  may  be  written 


(75) 
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defining 

CSL]  s /Cbnl]  T (Cs]  - [c]  [aeC]  ‘ [bnl]  dV 

V°  (77) 

and 

{pt}  * / Cbl]  T ({S}  ' [c]  {AE>C}  " {aeT})  dV 

V°  (78) 

yields  the  equations  of  motion  given  by  equation  (60). 

It  should  be  noted  here  that  an  analysis  similar  to  that  performed  In 
the  section  dealing  with  the  general  constitutive  law  will  show  that  the 
constitutive  tensor  for  this  loading  case  Is  Identical  to  that  obtained  In 
equation  (21). 

CONCLUSION 

1 

It  has  been  shown  by  Halsler  that  for  cyclic  load  histories  neither 
Isotropic  nor  kinematic  hardening  models  accurately  predict  the  Bauschlnger 
effect  for  many  materials.  It  has  been  shown  further  that  the  combined 
Isotropic-kinematic  hardening  rule  more  properly  depicts  material  response 
when  Isothermal  conditions  are  encountered.  A material  model  which  accounts 
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for  the  Bauschlnger  effect  and  Includes  both  creep  and  thermal  response  has 
been  proposed  In  this  paper.  It  Is  believed  by  these  authors  that  this 
model  will  more  accurately  predict  material  response  than  a kinematic  hard- 
ening rule  when  thermal  loads  are  present. 

The  proposed  model  presented  here  Is  currently  being  studied  In  the 
finite  element  code  AGGIE  I.  Results  of  this  study  are  forthcoming. 
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